'''
Created on Nov 30, 2010

@author: oabalbin
'''

#import os
import subprocess
#import pysam


def intervalsTorealign(ref_genome, indexed_bam_file, use_mem, num_cores, path_to_gatk, indeldb_file):
    '''
    time java -Xmx16g -jar /exds/sw/bioinfo/gatk/GenomeAnalysisTK-1.0.4705/GenomeAnalysisTK.jar 
    -T RealignerTargetCreator 
    -R /exds/projects/alignment_indexes/gatk/hg19/hg19.fa 
    -o /exds/users/oabalbin/projects/snps/exomes/aM18/test/realignment_output.intervals 
    -I /exds/users/oabalbin/projects/snps/exomes/aM18/test/s_3_12_sequence.hg19.aln2.rmdup.sorted.bam 
    -B:indels,VCF /exds/projects/alignment_indexes/gatk/hg19/dbsnp132_00-All_processed5.processed.vcf
    '''
    
    gatk_command=path_to_gatk+'GenomeAnalysisTK.jar'
    outfile=indexed_bam_file.replace('.bam','.realign.intervals')
    # '-nt', str(num_cores), not supported yet
    args = ['java','-Xmx'+str(use_mem)+'g', '-jar',gatk_command,'-T','RealignerTargetCreator',
            '-R', ref_genome, '-o', outfile, '-I', indexed_bam_file, 
            '-B:indels,VCF', indeldb_file]
            
    retcode = subprocess.call(args)    

    return outfile


def realign_bam_file(ref_genome, indexed_bam_file, realign_sites_file, indeldb_file, use_mem, num_cores, temp_dir,path_to_gatk):
    '''
    java -Xmx16g -Djava.io.tmpdir=/exds/users/oabalbin/projects/snps/exomes/aM18/test/temp 
    -jar /exds/sw/bioinfo/gatk/GenomeAnalysisTK-1.0.4705/GenomeAnalysisTK.jar -T IndelRealigner 
    -R /exds/projects/alignment_indexes/gatk/hg19/hg19.fa 
    -I /exds/users/oabalbin/projects/snps/exomes/aM18/test/s_3_12_sequence.hg19.aln2.rmdup.sorted.bam 
    -targetIntervals /exds/users/oabalbin/projects/snps/exomes/aM18/test/realignment_output.intervals 
    -o /exds/users/oabalbin/projects/snps/exomes/aM18/test/s_3_12_sequence.hg19.aln2.rmdup.sorted.realigned.bam 
    -B:indels,VCF /exds/projects/alignment_indexes/gatk/hg19/dbsnp132_00-All_processed5.processed.vcf -compress 0
    
    -B:indels,VCF /path/to/indel_calls.vcf \
    -D /path/to/dbsnp.rod \
    -knownsOnly \
    -LOD 0.4 \
    '''
    
    gatk_command=path_to_gatk+'GenomeAnalysisTK.jar'
    outfile=indexed_bam_file.replace('.bam','.realigned.bam')
    use_temp_dir='-Djava.io.tmpdir='+temp_dir
    #'-nt', str(num_cores), not supported yet
    args = ['java','-Xmx'+str(use_mem)+'g', use_temp_dir,'-jar',gatk_command,'-T','IndelRealigner',
            '-R', ref_genome, '-I', indexed_bam_file, '--targetIntervals', realign_sites_file,
            '-o', outfile,'-B:indels,VCF', indeldb_file, '-compress', '0']
            
    retcode = subprocess.call(args)
    
    return outfile


def intvTorealign_only_knownSites(ref_genome, use_mem, num_cores, path_to_intervals, path_to_gatk, indeldb_file):
    '''
    time java -Xmx16g -jar /exds/sw/bioinfo/gatk/GenomeAnalysisTK-1.0.4705/GenomeAnalysisTK.jar 
    -T RealignerTargetCreator 
    -R /exds/projects/alignment_indexes/gatk/hg19/hg19.fa 
    -o /exds/users/oabalbin/projects/snps/exomes/aM18/test/realignment_output.intervals 
    -I /exds/users/oabalbin/projects/snps/exomes/aM18/test/s_3_12_sequence.hg19.aln2.rmdup.sorted.bam 
    -B:indels,VCF /exds/projects/alignment_indexes/gatk/hg19/dbsnp132_00-All_processed5.processed.vcf
    '''
    
    gatk_command=path_to_gatk+'GenomeAnalysisTK.jar'
    outfile=path_to_intervals+'knownSites_realign.intervals'
    # '-nt', str(num_cores), not supported yet
    args = ['java','-Xmx'+str(use_mem)+'g', '-jar',gatk_command,'-T','RealignerTargetCreator',
            '-R', ref_genome, '-o', outfile, 
            '-B:indels,VCF', indeldb_file]
            
    retcode = subprocess.call(args)
    
    return outfile



def realign_bam_only_knownSites(ref_genome, indexed_bam_file, realign_sites_file, indeldb_file, use_mem, num_cores, temp_dir,path_to_gatk):
    '''
    java -Xmx16g -Djava.io.tmpdir=/exds/users/oabalbin/projects/snps/exomes/aM18/test/temp 
    -jar /exds/sw/bioinfo/gatk/GenomeAnalysisTK-1.0.4705/GenomeAnalysisTK.jar -T IndelRealigner 
    -R /exds/projects/alignment_indexes/gatk/hg19/hg19.fa 
    -I /exds/users/oabalbin/projects/snps/exomes/aM18/test/s_3_12_sequence.hg19.aln2.rmdup.sorted.bam 
    -targetIntervals /exds/users/oabalbin/projects/snps/exomes/aM18/test/realignment_output.intervals 
    -o /exds/users/oabalbin/projects/snps/exomes/aM18/test/s_3_12_sequence.hg19.aln2.rmdup.sorted.realigned.bam 
    -B:indels,VCF /exds/projects/alignment_indexes/gatk/hg19/dbsnp132_00-All_processed5.processed.vcf -compress 0
    
    -B:indels,VCF /path/to/indel_calls.vcf \
    -D /path/to/dbsnp.rod \
    -knownsOnly \
    -LOD 0.4 \
    '''
    
    gatk_command=path_to_gatk+'GenomeAnalysisTK.jar'
    outfile=indexed_bam_file.replace('.bam','.realigned.bam')
    use_temp_dir='-Djava.io.tmpdir='+temp_dir
    cleaning_threshold=0.4 # used in gatk guidelines
    #'-nt', str(num_cores), not supported yet
    
    args = ['java','-Xmx'+str(use_mem)+'g', use_temp_dir,'-jar',gatk_command,'-T','IndelRealigner',
            '-R', ref_genome, '-I', indexed_bam_file, '--targetIntervals', realign_sites_file,
            '-o', outfile,'-B:indels,VCF', indeldb_file, '-knownsOnly', '-LOD', str(cleaning_threshold),
            '-compress', '0']
            
    retcode = subprocess.call(args)
    
    return outfile



'''

for each lane.bam
    realigned.bam <- realign(lane.bam) [at only known sites]
    dedup.bam <- MarkDuplicate(realigned.bam)
    recal.bam <- recal(dedup.bam)

for each sample
    recals.bam <- merged lane-level recal.bam's for sample
    dedup.bam <- MarkDuplicates(recals.bam)


for each lane.bam
    realigned.bam <- realign(lane.bam) [at only known sites]
    dedup.bam <- MarkDuplicate(realigned.bam)
    recal.bam <- recal(dedup.bam)

for each sample
    recals.bam <- merged lane-level recal.bam's for sample
    dedup.bam <- MarkDuplicates(recals.bam)
    realigned.bam <- realign(dedup.bam) [with known sites if possible]

'''

























    